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ABSTRACT 


Response of composite laminates under in-plane compressive loading is of interest 
to the analysts and designers. Since they are generally thin, they are prone to insta- 
bility under inplane loads. One aspect which distinguishes composite materials from 
isotropic materials is the shear strength of the materials used.In the case of isotropic 
materials elastic modulus(E) is generally about 2. 5-2.8 times shear modulus(G), for 
composite materials since E is independent of G it can have any value.Hence the 
transverse shear effect becomes important for predicting the response of even moder- 
ately thick composite laminates.For very thick laminates neglecting transverse shear 
will lead to completely erroneous results.A number of different theories have been 
suggested by different investigators which take into account of this transverse shear 
effect. 

In the present study an attempt has been made to take into account the effect of 
transverse shear for the stability analysis of moderately thick/very thick compos- 
ite laminated plates under in-plane compressive and shear loading using a “Simple 
Higher Order Shear Deformation Theory” based on four unknown displace- 
ment functions instead of five which is commonly used for most of the other higher 
order theories.The finite element method is employed to study the initial buckling 
load of square/rectangular laminates.A continuous shear flexible finite ele- 
ment based on the proposed HSDT is developed using the Hermite cubic rect- 
angular element. The plate geometry has been modeled using linear Lagrange 
rectangular element and hence a Subparametric formulation has been adopted. 

The analytically obtained results using the present formulation have been compared 
with the results obtained employing different shear diformation theories.The change 
in initial buckling response of thick rectangular laminated plates with respect to 
the orientation angle has been studied.The interaction curves (between Nx and Ny 
and between Nx and Nxy for different aspect ratios) are studied in details. Several 
new results have been presented for different parameters like, aspect ratio, fibre 
orientation angle, a/h ratio, etc. 
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Chapter 1 


INTRODUCTION 


1.1 COMPOSITE MATERIALS 

Composites can be broadly defined at a macroscopic scale as a combination of two 
or more materials in which the materials remain in separate phases and they donot 
interact chemically. Common materials which can be included in this broad definition 
include reinforced concrete, fibre reinforced plastic, plywood, etc.Among the two 
phases one is a binder or matrix and the other is a reinforcement.The composite 
material discussed in the present work belongs to this category. 

A combination of dissimilar materials can have unique and very advantageous prop- 
erties if the materials have appropriate characteristics.lt may result in better prop- 
erties than either of the materials alone.For instance, the FRP pipe will be stronger 
aginst other traditional pipe materials such as all plastic PVC pipe.By forming a 
composite material one can improve a number of structural and behavioural proper- 
ties such as specific strength, specific stiffness, fatigue life, hygro-thermal behaviour, 
corrossion resistance, wear resistance, etc. 

Broadly, composites can be classifed into two types: 

a) Particle-reinforced composites 

b) Fibre-reinforced composites. 

In the present investigation laminated continuous fibre reinforced compos- 
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ite is considered.Some of the commonly used fibre materials are E-glass, S-glass, 
kevlar, graphite, boron, silica, tungsten, etc.Fibres are used to provide the nec- 
essary strength and stiffness to the structure.Matrix materials (Polymers) can be 
classified into two types depending on their structure and thermal behaviour.These 
are namely; 

1) Thermoplastic Polymers - Soften and melt on heating.Polyethylene, polystyrene, 
nylon, polyacetals belong to this group. 

2) Thermosetting Polymers - Do not soften but decompose on heating. Once 
solidified by a crosslinking (curing) process, they cannot be reshaped. Common ex- 
amples of this kind include epoxides, polyesters, phenolics, silicone, etc. 

The main advantages of polymeric matrix are low cost, easy processibility, good 
chemical resistance and low specific gravity. 

Composite offers unique opportunities in design.Apart from being lighter it is a good 
substitute for conventional metallic materials.The structural performance offered 
by it is far versatile than can be realised with conventional materials.An in-depth 
understanding of the principles governing their structural behaviour can help us 
explore this versatility to maximum possible extent. 

1.2 MOTIVATION 

As mentioned above the key factor in utilising the underlying strength and unique- 
ness of composite materials is the proper understanding of its structural response 
under different working load conditions.One structural component which is being 
extensively used in aerospace, automobile and other applications is laminated com- 
posite plates.Earlier, classical lamination theory which is an extension of the classi- 
cal plate theory was used to analyse laminated composite plate.Unfortunately this 
theory though simple, completely ignores the transverse shear effect which has sig- 
nificant effect on the structural behaviour of laminated plate.The reason behind this 
is the very high ratio of effective elastic moduli along the fibre direction to effective 
shear moduli which can be achieved in case of composite materials as the shear 
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moduli can be varied independent of elastic moduli though the same is not the case 
for isotropic materials.This when coupled with thickness effect makes classical lami- 
nate theory unusable for the analysis of even moderately thick composite laminated 
plates. 

One behavioural aspect of composite laminated plates which is vital for its safe 
design is the initiation of buckling under the action of in-plane compressive and shear 
loads.The load at which an initially fiat configuration of the plate under the action of 
in-plane forces suddenly changes to a bent configuration is known as initial buckling 
load.As the structure undergoes form failure at this point it is imperative that a 
considerable amount of attention must be paid to find its initial buckling load.The 
exact determination of this load governs the safe designing of plate structures. 

So, for a better estimation of the initial buckling of composite laminated plates, many 
theories have been proposed over the years which includes the effect of transverse 
shear deformation.The first order shear deformaion theory (FSDT) includes shear 
correction factor in their formulation, the exact evaluation of which is very difficult 
as the shear correction factor depends on the geometry and Poisson ratio.The higher 
order theories (HSDT) which take into account the parabolic variation of transverse 
shear stresses, normally involves five unknown displacement functions and requires 
(7° continuous element for a finite element analysis.This requires numerical strategies 
such as reduced/selective integration to avoid shear locking effect (see [1]).A higher 
order theory proposed by Lim et al. [2] involves only four unknown displacement 
functions which allows continuous plate bending element and therefore doesnot 
necessitate numerical strategies to avoid shear locking eflFect.Also as this theory in- 
volves lesser number of unknowns it is expected to be computationally more effective 
than the other theories which involve more (five) unknowns. 


1.3 LITERATURE SURVEY 

Over three decades analysis of laminated composite plates is under intensive inves- 
tigation.Some of the studies ([3] and [4]) have shown that the transverse shear effect 
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is quite significant in the layered composite plates due to high ratio of inplane elastic 
modulus to transverse shear modulus which made classical laminate theory unusable 
for analysis. To overcome this difficulty a modified bending theory was formulated 
which accounts for the transverse shear deformation in isotropic plates. This plate 
bending theory is known as Reissner/Mindlin plate bending theory ([5] and [6]).This 
theory was applied for anisotropic plates by Whitney and Pagano (see [3] and [4]). In 
this theory the transverse shear strains are constant throughout the thickness of the 
plate. A shear correction factor has to be introduced to account for the discrepency 
between the theory and the actual behaviour.Though this theory gives better re- 
sults than classical laminate theory, however, gives erroneous results when the plate 
thickness to side ratio increases.Thus to account for a better representation of shear 
distribution across the thickness, Reddy [7] proposed a higher order theory.Reddy 
and his co-researchers studied the initial buckling of laminated composite plates us- 
ing this higher order shear deformation theory ([8] and [9]). A simplified higher order 
theory proposed by Lim et al. [2] involves only four unknowns instead of five for 
Reddy [7]. The simplified theory allows to use continuous plate bending element 
which is free from any shear locking effect.Gajbir Singh [10] used this model for 
the response of laminated composite beams and plates.His work was related to non- 
linear bending, vibration and thermal postbuckling of composite plates.To the best 
of the author’s knowledge no results are available for the initial buckling analysis 
under uniaxial and biaxial compressive and shear inplane loading in the published 
literature using this simple higher order shear deformation theory. 


1.4 PRESENT INVESTIGATION 

In the present investigation the following studies have been made: 

1. Initial buckling behaviour of laminated composite plate under in-plane com- 
pressive loading using the simple higher order shear deformation theory as mentioned 
in section 1.3 and compare the results of other investigators using other higher order 
theories. 

2. The effect of various parameters like, plate aspect ratio, length to thickness 
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ratio and lamination angle on tlie initial buckling load of laminates. 

3. Development of interaction curves for biaxial loading {N^ and Ny, Nx and 
+Nxy, Nx and —Nxy) of symmetric laminate for different aspect ratio, length to 
thickness ratio and lamination angles. This will help in the design of laminates 
subjected to combined loads. 


1.5 LAYOUT OF THESIS 

The organisation of the thesis is given as: 

1. The first chapter is devoted to the motivation for present study, the lit- 
erature review and the goals pertinent to the present study and the layout of the 
thesis. 

2. The second chapter is devoted to the discussion of the plate model and 
laminate constitutive equations. 

3. The third chapter deals "with the finite element formulation of the problem 
using energy principle. 

4. In the fourth chapter, fi.rst the validation of the present model is presented 
and afterwards various numerical results are furnished to detail the present investi- 
gation and inferences on them are drawn. 

5. Concluding remarks and observations from this study along with some sug- 
gestions for further extensions of the work are given in Chapter five. 
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Chapter 2 


HIGHER ORDER PLATE 
THEORY 

2.1 INTRODUCTION 

In this chapter, details of a plate theory proposed by Lim et al. [2] for isotropic 
plate is discussed and extended to laminated composite plate.The displacement field 
includes classical plate theory and first order shear deformation theory as its subset 
and accounts for the parabolic variation of the transverse shear strains which also 
accounts for the surface boundary condition of zero transverse shear stresses (and 
hence shear strains) at the top end bottom surfaces of the plate. 

The displacement field presented by Reddy [7] is modified by reducing the number 
of variables to four (five, in case of Reddy’s theory) .This simplification is based on 
the assumption that in-plane rotation tensor is constant through the thickness and 
the transverse displacement can be seperated into two components in such a way 
that the transverse shear strains are functions of only one of these components. 
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Figure 2.1: Laminated plate with orientation of global co-ordinates 

2.2 FORMULATION OF DISPLACEMENT 
FIELD 

The displacement field is expressed as: 

{V (a;, y,z)} = [ u{x, y, z), vix, y, z), w{x, y, z) f . 

where u{x,y,z), v(x,y,z) and w(x,y,z) are: 

u{x, y, z) = uo(x, y) - zwb,Jx, y) + z^(j>x{x, y) + y)] 

v(x,y,z) = vo{x,y) — zwb,y{x,y) + z^(j)y{x,y) + z^ipy{x,yy, (2-1) 

w{x, y, z) = Wb{x, y) 4- 'Ws{x, y) 

The transvBrsG displa-CGiiiGiit compoiiGiit is such that its dGrivativGS arc numeri- 
cally equal to the rotation of the cross-section (i.e. (j) = -sjwi) and Ws is the dis- 
placement due to the effect of transverse shear deformation of the cross-section. The 
assumption that W}, and Ws are functions of x and y only is justified, since the 
transverse normal stress is of the order of (h / times the in-plane normal stress. 


7 




The von-Karman type non-linear strain-displacement relationships for the higher 
order theory under consideration can be written as follows: 

^xx = Uo,x + I {Wb,^ + Ws,xy - Z -f 

% = ^0,J/ + I {Wb,y -1- Ws,yy - 2 : Wb,yy T (j)y,y + , y, 

^ZZ ~ 0) 

Ixy = Uo,y + Vo,^ + {Wb,^ + Ws,^)(Wb,y + Ws,y) - 2z Wb,^y (2-2) 

-\-Z (pxyy) Z fpxiy)] 

Ixz = Ws,^ + 2 Z (f)^ + 3 z‘^ Ipx] 
lyz = '^a,y + 2 Z (j)y + Z Z^ tpy 

where comma (,) denotes the partial derivatives. 



Figure 2.2: Representation of arbitrary three dimensional domain 

The surface boundary condition that the transverse shear stresses vanish on the 
plate’s top and bottom faces (see Fig. 2.2) is equivalent to the requirement that the 
corresponding strains be zero on these surfaces, i.e. 

Txz(^5 2/) ^ 2) ~ 'Yyzi^j ^2) ^ 

On introduction of the conditions as given above in the expressions for transverse 
shear strains (from Eq. 2.2), the following relations are obtained. 

(f)x = <l)y = 0 and 

'^x — ~3h^ Ws,x ] ' 02 / ~ ~~3h^ '^^>y 


(2.3) 



The displacement field of Eq. (2.1) is modified by setting 4>^ and (j)y to be zero 
according to conditions of Eq. (2.3). The resulting displacement field is written 
below: 

u{x, y, z) = uo{x, y)-z Wb,^{x, y) - ^ Ws,x{x, y); 
v{x, y, z) = vo{x, y)-z Wb,y{x, y) - ^ Ws,y{x, y); 
w[x, y, z) = Wb{x, y) + w^ix, y) 

In Eq. (2.4) u, v and w are the displacements along x, y and z directions respectively- 
uo,vo and w {wb and Ws) are the mid-plane displacements.Thus, the generalised 
displacement vector {5} of the mid-surface contains four degrees of freedom (DOF) 
and is given by: 


{ 5 } = {uo,vo,Wb,Wsy. 

The corresponding strain-displacement relationships are: 

^xx = Uo,x + I - 2: Wb,^x - Ip fj^s,xx'> 

tyy = VQ,y -I- \ {Wb,y + Ws,y)^ “ ^ Wb,yy “ Ws,yy] 

^zz = 0; 

7xy = uo,y + no,x + {Wb,x + 'Ws,x){wb,y + Ws,y) -2z Wb,xy 

- nij 

3/i2 ^s^xy'i 

'Jxz — (1 ~ "p) resjx) 

'Yyz = (1 - p) Ws,y 

It can be noted that these kinematics can be obtained by substituting tpx = —'^b,x^ 

- Wb,y and w = Wb + Ws in Reddy’s [7] higher order theory involving five unknowns 
{Uo,Vo,Wo,'lpx,i^y)- 

It can also be shown that the present formulation degenerates to first order shear 
deformation theory by suitably dropping higher order terms and incorporating ap- 
propriate shear correction factors (for isotropic plate whose value is 5 / 6 ) -The clas- 
sical laminate theory (CLT) can be derived from the present theory by equating n’* 
to zero. 
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The strain-displacement relationships given in Eq. (2.5) can be rewritten as: 


^xx = £0x + Z Kqx + Klx, 

€.yy = eojr + 2 Koy -+• Kiy‘, 

^zz ~ 0 ; 

Ixy = 'JOxy + Z KQxy + Z^ ^Ixy', 

Ixz = (1 ~ 

lyz = (1 — %■) Wsiy 


where: 


eox = no,, + i {Wb,^ -I- 

^Oy “I" 2 {^b^y T n^53y) 3 

iQxy = Uo,y + n03x + (^63 x + Ws,a:){Wb,y + Ws,y)l 

^Ox '^btxxi ^Oy — ^bjyy] ^Oxy ~ ‘^'^bixy'i 

4,4 8 

f<'lx — ““"sp’ '^S’txx') ^ly ^ ~3p '^sjyy] ^Ixy = "^sp '^Sixy 

2.3 LAMINATE CONSTITUTIVE 
EQUATIONS 


A lamina (considered a unidirectional composite) is characterised by having all fibres 
(either a single ply or multiple plies) oriented in the same direction.This model 
allows one to treat the lamina as an orthotropic material whose material symmetry 
planes are parallel and transverse to the fibre direction.The model used to represent a 
lamina consists of unidirectional fibres per layer which is further assumed as perfectly 
straight and uniformly oriented within the lamina. 

In the following Fig. (2.3) the material co-ordinate axes L and T are defined parallel 
and perpendicular to the fibre direction respectively, while global co-ordinates are 
X, Y and Z. The material co-ordinate axis T' is further assumed to be along Z 
axis.The angle between global co-ordinate axis X and fibre direction, material co- 
ordinate axis L, is 9 and is known as the orientation angle. As a sign convention 
anticlockwise angle from +ve X direction is taken as positive. 
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Figure 2.3; Co-rdinate axes for the leimina 
Generalized Hooke’s law 

In the derivation of lamina constitutive equations, following two assumptions are 
made. 

1. The lamina is a continuum. 

2. It behaves as a linearly elastic material. 

Further, at the micro-level the following assumptions are made about the material. 

1. Perfect bonding between fibres and matrix exists. 

2. Fibres are parallel and uniformly distributed throughout. 

3. The matrix is free of voids or micro cracks and initially in a stress free state. 

4. Both fibres and matrix are isotropic and obey Hooke’s law. 

Stress-strain relations for the lamina in the material coordinate axes, whose fibres 
are oriented at an angle 6 with reference to the X axis is given as (see [11], [12]): 

{ a, }^ = [Qy]fc { e, (2-7) 

where is the vector of stress components, [Qij] is the 2D- stiffness matrix 
with respect to material co-ordinate axes, and {^j} ^-re the engineering strain com- 
ponents, for the lamina. 

In expanded form, the above relation can be written as: 
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where: 

Q22 = (i_„ 5 ^^)i Om = Glt; 
Q 44 ~ ^lt' 'i Q 55 ~ G'pj^ 


The stresses and strains in the X, Y and Z directions are obtained by transformation 
of the relations given in equation (2.7). The transformed stress strain relations for 
the lamina are given as: 

f 


‘xy 

'^XZ 

'^yz 


> = 


On 

O 12 

016 

0 

0 


f \ 

O 12 

O 22 

026 

0 

0 


Cy 

O 16 

O 26 

Qm 

0 

0 

< 

7 xy ' 

0 

0 

0 

O 44 

045 


'Jxz 

0 

0 

0 

O 45 

055 

k 



( 2 . 8 ) 


where: 

011 = OllC^ + (2 Qi2 + 4Q66)®^C^ + 022^^) 

012 = <512(5^ + + (Oil + O22 ~ 4066 )^^^^) 

O16 = (On - O12 - 2 Qm)sc^ + (O12 - O22 + 2 Q 66 )s^c; 
O22 — Oii®^ T (2O12 + 4066 ) 5 ^*'^ "h O22C ; 

O26 = (On ~ O12 — 2Q66)s^c + (O12 — O22 + 2 Q 66 )sc^; 
Oee ~ (On ~ 2 ( 3 i 2 + O22 ~ 2066 )s^c + Q 6 s{c + s ); 

0 44 = 044^^ + Oss^^j 

045 “ (O44 — Oss)^^; 

O55 = OssC^ + 044 ^^ 
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where, [Q^j] is the transformed stiffness matrix , c = Cos{9) and s = Sin(9) with 9 
denoting the angle of orientation of lamina. 

Using the above lamina constitutive equations and integrating the stresses over the 
laminate thickness, the stress resultants in terms of inplane forces, moments, shear 
forces and higher order stress couples per unit length, for a laminate with NL 
laminae are, 

' {Ni} 1 r [A,,] [B,,] [E,,] [0] 1 f {60,} ' 

^ {Mi} 1 [B,,] [A,] [F.,] [0] I {«o,} I ^2 9) 

{Pi} I [A,] [A,] [A,] [0] I {«i,} I 

^ {Ql} ^ _ [0] [0] [0] [Gim] _ ^ {7m} ^ 

(i , j = 1 , 2 and 6) 

(1 , m = 4 , 5) 


(Detailed expansion of Eq. (2.9) is given in Appendix A.) 

Where {eo} are the mid-surface strains, {/cq} are the mid-surface bending and twist- 
ing curvatures, {ki} are the higher order terms and {7} are the transverse shear 
slope and 


r^k-i 

{Ni}’^ = {Nx, Ny, N^cy) = 

fe=l '’’‘k 

{M!i} = My^ Mxy) — ^ } / iPxi^yi'^xy) ^ 

k=l '’’‘k 


m'^ = {Px,Py,Px 


NL » 

,) = E/ 

L ._1 J Z\ 


(c^a;, O'j^, L'xy') Z dz^ 


fc=l •'^k 


( 2 . 10 ) 


( 2 . 11 ) 


( 2 . 12 ) 


iA r^k-i. 4^2 

{Ql}^ = {Qx, Qy)=J2 (2.13) 

i:=l •'^k 


Here, z^ are the z coordinates corresponding to the laminae interfaces as shown in 
Fig 2.4. 


Thus, with the assumed displacement model, the various stiffness matrices derived 
are: 

[Ay] = Extensional; [Aj] = Bending-extensional coupling; 
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Figure 2.4: Geometry of multilayered laminate 

[D^j] = Bending; [Gim] = Transverse shear stiffness and 
[Eij], [Fij], [Htj] — Higher order stiffnesses. 


The set of stiffness matrices [Aij\, [By], [By], [By], [Fi^], [By] and [Gjm] are used in 
forming overall stiffness matrix [Dr] for the laminate. 
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Chapter 3 


FINITE ELEMENT 
FORMULATION 


For the finite element approximation of the problem under investigation rectangular 
four noded elements are used, along with linear Lagrange interpolation function to 
model the geometry of the plate. A continuous shear flexible element based on the 
presented higher order theory is developed using the Hermite interpolation formulae 
as indicated by Bogner et al [13]. A typical mesh, generated over the rectangular 
plate domain is shown in Fig. 3.1. The node (both local and global) and element 
numbering scheme is also shown in the figure. In what follows the finite element 
formulation is presented. 


3.1 DEFINITIONS 


Let {y} be the displacement vector defined as: 


r \ 

u 


{V} = { V 


W 

\ / 


(3.1) 
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B 



Element node numbering 
of i-th Element 

5 and Ti are the natural co-ordinate axes of 
master element 

X and Y are the local co-ordinate axes of 
physical element 


Figure 3.1: Typical meshing of Rectzingular plate domain 


16 



The stress and strain vectors corresponding to {F} be {a} and {e} can be defined 



From the generalized Hooke’s law which relates stress components to the respective 
strain components in global coordinate system: 

{a} = [Q] {6} (3.3) 

The material stiffness matrix [Q] for orthotropic material is given in Eq. (2.8) for 
each lamina. 

We are concerned with initial buckling of the laminate.Hence, the non-linear terms 
in the strain-displacement relationships as described in Eq. (2.5) are neglected for 
our study because the transverse deflection (both Wb and Wg) is assumed to be very 
small. 

The linear strain-displacement relationships can be written as: 

Cxx “ ^ '^h-ixx ~~ '^s^xx^ 

4tZ^ 

^yy — ~ ^ ^f>’yy ~ 

= 0 ; 

Try ~ “b ^0)a: ^ ^ '^bjxy (^■^) 

— Mi -in 

3 / i 2 '^sixyi 

Ixz “ (f ~JjT ) 
lyz — (1 ~ ) ^Siy 

The above relations is expressed in matrix form as: 

{e} = [D] {F} (3.5) 
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where [£>] is a differential operator in terms of global coordinates, defined as 


[D] = 


dx 

0 

0 

0 

dy 

0 


d_ 

0 

dy 

dx 

dz 

0 

d_ 

dx 

0 

d 


dz 

dy 


(3.6) 


The components of displacement can be expressed in terms of the four unknowns, 
which can be written as: 


{5y = {uo,Vo,Wi,Ws} 


For any element the field variables in terms of shape functions and nodal variables 
can be written as: 

4 4 4 

Wo (a;, y) = woj + 53 -^*+4 + 53 

l—X 2=1 2=1 

4 4 4 

Wo(a;, y) = 53 ^0* + 53 ^*+4 + 53 -^+8 

2=1 2'=1 2=1 

4 4 4 

y) = 53 + 13 -^*+4 + 53 ^^'+8 ^ b,y + 

2=1 2=1 2=1 

4 

“^ 2+12 

2=1 

4 4 4 

Ws{x, y) = ^ iV' ii;52 + X!' ^^2+4 yJs,xz + -^ 2+8 ^ 5 , 2 ,^+ 

2 *=! 2=1 2*=1 

4 

53 ^ i +12 
2=1 



In the above expression uo„ uo,y., voi, vq,^^, vo,y., wi,i, Wb,^^, Wb,y^, “^bixyii 
Ws,xi^ '’^s,y^, '^s,xy^ (* = are the fourteen degrees of freedom per node, N'^s are 
the element shape functions in local co-ordinates (in terms of X and Y as shown in 
Fig. 3.2) .In natural co-ordinates (in terms of ^ and rj as shown in Fig. 3.2) N[s can 
be expressed as N^s (see [l]).The details of the element shape functions (both iV( 
and iV() are given in Appendix B. 
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Y 



(N,',N^N^N,'3) ® O^N^N,'oN,'4 


Four noded rectangular element in 
local co-ordinate system 





Four noded square element in 
natural co-ordinate system 


Figure 3.2; Four noded rectangular element 

Substituting Eq. (3.7) in the strain-displacement relations given in Eqn. (3.4) the 
following relations can be established: 



where: 

{ 5 }^ = { Uqi, liOjj/j) ‘^Oij '^biyp “^b-ixy^i '^arixii 

'Ws,xy^}> (* = l.->4) 

and: 
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^Ox ^Oy — IfOxy — + ^Ojxj 

f^Ox '^b'ixx') ^Oy ^^^yy'i ^Qxy ~ ^'^b-ixyi 
'^Ix Ws,xx> '^ly ~ ~'^ '^Sjyy', ^^Ixy = — Wg^^y 

The details of [Bq] are given in Appendix C. 

3.2 FINITE ELEMENT FORMULATION USING 
ENERGY PRINCIPLE 

The total potential for the plate is given by 

nW=Eir'W, (3.9) 

e=l 

where tt® is the total potential of the non-intersecting (but adjacent) sub-domains 
e which are part of the domain (N sub-domains are considered here). The total 
potential can be expressed in terms of internal strain energy external work 
done due to transverse loads and work done due to inplane forces T^®\ as 
follows; 

TT® ( 5 ) = f/^®) - w^®) - r(®^ 

Strain energy of the laminate can expressed as follows: 

= \ f {e}^ [a] dF(®) 

2 JvM 

Work done by the applied external transverse load is, 

= /a+w » w f- 

where /+ is the transverse load on the top face and f~ is the transverse load on 
the top face R~ and w = wi) + Wg. 

For initial buckling problem we can assume that the plate, which is stressed by 
forces acting in the middle plane of the plate (i.e. Ny and N^y), undergoes 
some small lateral bending consistent with the given boundary conditions and in- 
plane forces remain constant during bending.Such limited bending can be produced 


(3.10) 

(3.11) 

(3-12) 
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Figure 3.3; The orientation of negative in-plane forces 

without stretching of the middle plane.Hence the work done by the applied in-plane 
loads (considered negative in sense) can be expressed as follows: 

= ]■ [ [A(j, w,l + w,l + 2 Na;y w,^w,y] dA (3.13) 

Where the directions of the inplane forces are as shown in Fig 3.3 . It can also be 
noted here that in case of initial buckling problem the applied external transverse 
load is kept zero (i.e. = 0). 

The exact solution {Vex} to this problem is the minimizer of the total potential H. 
This can be obtained as: 

5(1) n = 0 (3.14) 

which is also the ViHual Work Formulation of the problem in terms of the com- 
ponents of 5 (and hence 5). This leads to the generalised finite element formulation 
which can be expressed as: 

[K] {5} = {F} (3.15) 

Where [K] is the global stiffness matrix which is composed of two components (i.e. 
global elastic stiffness matrix [K] and global geometric stiffiiess matrix [ifg]) and 
[F] is the global load vector which is zero in our case as W(®) = 0. 
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Hence Eqn. (3.15) can be rewritten as: 


[^1 - IK,] { 5 } = 0 


(3.16) 


In the following sections, we are going to employ the virtual work formulation to 
derive the finite element equations governing the initial buckling problem. 


3.3 COMPUTATION OF ELEMENT ELASTIC 
STIFFNESS MATRIX 

The elastic stiffness matrix corresponding to assumed deformation state of an ele- 
ment can be defined by expressing the internal strain energy in terms of unknown 
nodal displacements. In the formulation of unsymmetric laminates the membrane, 
the flexure, membrane-flexure coupling and shear strains contribute to strain energy. 

Prom Eqn. (2.9 - 2.13) and Eqn. 3.11 the internal starin energy of an element (given 
by area A®) can be expressed as: 




( Nx ^Qx d" ^Qy "f" iV py €()xy d“ -^x ^Ox d” ^Oy d" ^^xy ^Oxy d" Px ^Zxd“ 
Py l^ly d" Pxy l^lxy d" Qx ‘^s,x d" Qy '^s,y )dA 


(3.17) 


The Eq. (3.8) along with Eq. (2.9) can be used to express the internal strain energy 
as: 


= i / ({J'riVl[AlISol{i*})<M (3.18) 

^ Jaw 

The strain energy expression can be written in a concise form as: 

= i {5^Y [K^] {j®} (3.19) 

where [iiT®] is the elastic stiffness matrix and {5®} are the unknown displacements 
corresponding to the finite element solution in an element e. 
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[K^\ is given by; 

[^1 = Lw [So)’’ [D,l [B„l dA (3.20) 

Numerical integration is carried out to get the element elastic stiffness matrix. 

3.4 COMPUTATION OF ELEMENT GEOMET- 
RIC STIFFNESS MATRIX 

The total work done by in-plane applied loads in an element can be expressed as: 

r* = 5 / ( {«T [B/1 [JV] ft] OT ) (3.21) 

This can be written in a concise form as: 

= I [K^] {5®} (3.22) 

The details of [Bg] are given in Appendix C. 
whereas: 

A(e iV™ 

[AT] = 

Nxy Ny 

Assuming that forces N^, Ny and N^y are represented by certain expression with a 
common factor A, so that 

Nx — A Nx^ Ny = A Ny, and Nxy — A Nxy 

then: 

Nx Nxv 
[iV] = A / 

_ Nxy Ny 

or: 

[JV] = A IJV) 
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where: 


[N] = 


Nxy 

Nxy Ny 


[Kg] is the geometric stifiFness matrix and {J®} are the unknown displacements cor- 
responding to the finite element solution in an element e. 


[Kl] is given by: 

[K;] = A [ [Bgf [iV] [Bg] dA- 
Or [Kl\ = A [ifj] 


(3.23) 


Numerical integration is carried out to get the element geometric stifihess matrix. 


3.5 METHODOLOGY FOR FINDING INITIAL 
BUCKLING LOAD 


In the present analysis the goal is to find the initial buckling load for laminated com- 
posite plate under different conditions of loading (uniaxial or biaxial) with different 
boundary conditions, plate aspect ratios, plate thickness ratios, fibre orientation 
angles, etc. 

This is an eigenvalue type of problem which is described below. 

Prom Eqns. (3.16), (3.20) and (3.23) we get the final FEM equation as follows: 

[[IT] - A[ir, )]{?}= 0 (3.24) 

The characteristic equation of the above eigenvalue problem can be expressed as: 


[K] - \[Kg] 


= 0 


(3.25) 


Eqn. (3.25) is solved using NAG ROUTINES to obtain the eigenvalues A,- (i = 
1,M) where M is the dimension of the matrices [K] and [Kg] (and is equal to the 
total degrees of freedom in the analysis set of the problem) .The minimum value of 
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A corresponds to the scaled initial buckling load from which we can get the initial 
buckling load as follows: 

Nx ( or Ny or N^y ) 

^min Nx ( or Ny or Nxy ) 

The corresponding eigenvector (which is also found using NAG ROUTINES) gives 
the mode shape (not the exact displacements) corresponding to the initial buckling 
load. 

3.6 GEOMETRIC APPROXIMATION 

3.6.1 LINEAR MAPPING OF STRAIGHT EDGE RECT- 
ANGULAR ELEMENTS 

As mentioned earlier, the plate domain is divided into a number of straight edge 
rectangular elements. For the aforementioned geometry the geometric transforma- 
tion from physical element to the master element can be done by linear mapping 
which is discussed below. 


LINEAR MAPPING: 


(X4,Y4) 


(X3,Y3) 


(Xi.Yi) 


(X2,Y2) 



Master Element 


Physical Element 

Figure 3.4: Linear mapping in two dimensional domain 
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The transformation between {X,Y) co-ordinates and (^, 77 ) co-ordinates is ac- 
complished by the following co-ordinate transformation; 


J = 1 
4 




(3.26) 


j=i 

Where X® and are the global co-ordinates of the physical element nodes and 
7 /)® denote the finite element interpolation functions of the master element.In the 
present study bilinear Lagrange interpolation functions are employed. 


7 ^® are defined in the following way: 

= ^(1 - 0 (1 - ; 

= J (1 + 6 (1 - ??) ; 

^3 = '^(1 + 0 (1 + ^) ; 

k = j(i - 6 (1 + v ) 


(3.27) 


It is obvious from the above relations that: 

^ i [— Xi (1 — 77) + X2 (1 — 7?) + Xz (1 -I- 77) — X4 (1 + 77) ] ; 

= 1 [_ X: (1 - e) - ^2 (1 + e) + ^3 (1 + e) + ^4 (1 - 0 ] ; 

^ = i [- n (1 - ,) + n (1 - ») + n (1 + I)) - I4 (1 + I?) 1 ; 

= 1 [_ Ft (1 - 0 - n (1 + 0 + Fj (1 + f) + ii (1-01 

dr] 4 

(3.28) 


Transforming dxdy into d^dr]'. 

Let the differential area dxdy is formed through vectors dx and dy with magnitude 
dA and direction normal to the elemental area is k 


dA = dxdy = {dx x dy]-k 


(3.29) 


dxdy = {{^di''i + ^dri3)x{^d^iA^dr]3)]-k 

= \J\ de dr] 
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or the Jacobian matrix is given by 


[J] = 


dx dy 

dx dy 

L St? St? J 


or: 


\J\ = 


dx dy dy dx 


d^ dy d^ dr) 
Hence f (x, y) dx dy can be written as 


^ / {x, y)dxdy = f (^, 77) d^ dy 


(3.31) 


(3.32) 


(3.33) 


such that / (^, y) = \J\ f {x{^, y), y{^, y)). 

The condition which must to be satisfied for physically meaningful mapping of the 
physical element to the master element is that, \J\ > 0. 


Numerical Integration; 

To compute f (^, y) d^ dy Gaussian quadrature for two-dimentional integrals 
is used through which 

» M N 

/ / (^5 v ) dy « 

JAm j-x 


where, M and N denote the number of quadrature points in ^ and y directions 
respectively, (^t, y^) denote the Gauss points, Wi and Wj are the Gauss weights at 
quadrature point {^t,yj) and geometrically the Jacobian |J| is the ratio of an area 
element in the physical element to the corresponding area element in the master 
element e. 

Order of numerical integration is greater than or equal to (p + l)/2 where, p is the 
order of the integrand. 

In our problem we have chosen 5X5 Gauss quadrature points. 


3.6.2 SHAPE FUNCTION DERIVATIVES 

Shape functions are defined in the master element i.e. (^, y). Its derivatives 

with respect to global co-ordinates (X and Y) are related to their derivatives with 
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respect to the natural co-ordinates and rj) by: 


1-^ Kx 


( 3 . 35 ) 




y,l 2 x,^ y,^ 

y?r) 2 x,^ y,n 

^>v V’i V’n Vr 




y^^ri 



( 3 . 36 ) 
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Chapter 4 


RESULTS AND DISCUSSION 


4.1 INTRODUCTION 

The finite element formulation presented in Chapter 3, is used to study the initial 
buckling of laminated composite plates using a simple higher order shear deforma- 
tion theory.For the present study it is assumed that there is no coupling between 
the in-plane and out-of-plane quantities.If this is not true, then even for an infinites- 
imally small inplane compressive force the laminated composite plate will undergo 
flexure. For symmetric laminate this condition is automatically satisfied but for an- 
tisymmetric laminate there is always coupling between in-plane and out-of-plane 
quantities.In the case of antisymmetric laminate we have assumed that this cou- 
pling effect is negligibly small and was ignored to find the initial buckling load.This 
assumtion although seems strong it can be shown that if the number of plies is 
increased the coupling effect dies out rapidly and the solution approaches the or- 
thotropic solution (see [12]). 

In the present study the following assumptions are made. 

1. The thickness of each layer is constant. 

2. There is no delamination/slip between the layers. 

3. No micro-buckling takes place before the initial buckling. 

4. There is no fibre breakage or matrix cracking before the onset of initial buck- 
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ling. 

In the current investigation the following nomenclatures are assumed to denote the 
thickness of laminate. 

• a//z > 100 : Thin laminate . 

• 20 < a//i < 50 : Moderately thick laminate . 

• a//i < 10 ; Very thick laminate . 

This classification is tentative and actual definition depends on a variety of factors 
such as loadings, boundary conditions etc. 

4.2 BOUNDARY CONDITIONS 

The following types of geometric boundary conditions are considered in the analysis. 

• Simply-Supported (S) : Wb = Wg = 0 . 

• Clamped (C) : Wb = Ws = Wb,^ = = 0 ; 

where, n = x ox y depending on the side of the plate. 

To define the complete set of boundary conditions on all four sides of the plate the 
following scheme has been adopted which is explained through an example. 

A set of boundary conditions say, SCCS denotes the following: 

Wb = Ws = 0 at a: = 0 V y 

Wb = Ws = Wb,y = Wg,y = 0 aty = 0Va; 

Wb = Ws = Wb,^ = = 0 ata; = aVy 

Wb = Wa = 0 dXy = h^x 

4.3 MATERIAL PROPERTIES USED 

For the present investigation the following material properties are used. 
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Table 4.1: Material properties 


Material code 

EiiOPa) 

Mjl 

Ex 

Qml 

Et 

LT 

Ex 

^tt' 

Ex 

^LT 

Ml [14] 

130 

13 

0.5 

0.5 

0.2 

0.35 

M2 [9] 

200 

40 

06 

0.6 

0.5 

0.25 

M3 

100 

isotropic 

0.25 

M4 [15] 

153 

15.904 

0.616 

0.616 

0.336 

0.32 


4.4 CONVERGENCE STUDY 

In this section the following set of problems are considered for the convergence study. 


Table 4.2: Problem definitions for convergence study 


Problem 

number 

Material 

code 

Aspect 

(a/b) 

Length / thickness 

ratio(a//i) 

Boundary 

conditions 

laminae lay-up 

sequence 

PI 

Ml 

1.0 

10 

ssss 

[90/0], 

P2 

Ml 

1.0 

50 

ssss 

[90/0], 

P3 

Ml 

1.0 

100 

ssss 

[90/0], 

P4 

Ml 

1.0 

10 

cccc 

[90/0], 

P5 

Ml 

1.0 

50 

cccc 

[90/0], 

P6 

Ml 

1.0 

.i 

100 

cccc 

[90/0], 

P7 

Ml 

— 1 

1.0 

10 

ssss 

[90/0/90/0]2 

P8 

Ml 

1.0 

1 

50 

ssss 

[90/0/90/0]2 

P9 

Ml 

1.0 

1 

100 

ssss 

[90/0/90/0]2 

PIO 

Ml 

1.0 

10 

cccc 

[90/0/90/0]2 

Pll 

Ml 

1.0 

50 

cccc 

[90/0/90/0)2 

P12 

Ml 

1.0 

100 

cccc 

[90/0/90/0)2 


Convergence of the initial buckling load with mesh refinement (by decreasing the 
mesh size or in other words by increasing the number of meshes in the plate domain) 
for the problems defined in Table 4.2 is as shown in Table 4.3 . 
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Prom Table 4.3 it can be concluded that for a square plate, 10 x 10 mesh gives suffi- 
cient accurate results (for both simply supported and clamped boundary conditions) 
from engineering point of view. Hence for further studies of square plates, the plate 
domain is divided into 10 x 10 mesh. 

For rectangular plates, convergence study was done for each problem (with different 
aspect ratios) seperately, trying to keep the element aspect ratio almost equal to 
one.It has been observed that for a/h = 2 sufficiently converged result can be ob- 
tained for 10 X 5 mesh whereas for a/h = 3 it is 15 x 5, for both simply supported 
and clamped boundary conditions. 


Table 4.3: Convergence study 


Problem number 

Mesh size 

Nondimensionalized buckling load, 

(Nx = N^h^/Eth?) 

PI 

2x2 

7.4943 

4x4 

9.3565 

8x8 

9.4372 

10 X 10 

9.4441 

P2 

2x2 

9.9005 

4x4 

12.4832 

8x8 

12.5253 

lOx 10 

12.5259 

P3 

2x2 

10.0013 

4x4 

12.6150 

8x8 

12.6554 

10x10 

12.6554 

P4 

2x2 

12.6457 

4x4 

15.3898 

8x8 

15.9465 

lOx 10 

16.0501 
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Problem number 

Mesh size 

Nondimensionalized buckling load, Nx 
{Nx = Nx9/Eth?) 

P5 

2x2 

20.3343 

4x4 

26.1606 

8x8 

27.4415 

10 X 10 

27.7342 

P6 

2x2 

20.6987 

4x4 

26.7411 

8x8 

28.0754 

lOx 10 

28.3809 

P7 

2x2 

10.3057 

4x4 

10.9218 

8x8 

11.0400 

10 X 10 

11.0477 

P8 

2x2 

13.1137 

4x4 

13.4850 

8x8 

13.4996 

10 X 10 

13.4997 

P9 


13.2234 


13.5846 

8x8 

13.5957 

lOx 10 

13.5962 

PIO 

2x2 

18.9513 

4x4 

21.2392 

8x8 

21.7043 

lOx 10 

21.7440 

Pll 

2x2 

39.8804 

4x4 

42.6732 

8x8 

44.5310 


44.7049 
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Problem number 

Mesh size 

Nondimensionalized buckling load, 


2x2 

40.7152 

P12 

4x4 

43.8164 


8x8 

45.8574 


10 X 10 

46.0485 


4.5 VALIDATION STUDY 


The validation of the plate model for finite element formulation is done by comparing 
the results of the following set of problems with standard results published in the 
literature. 


Table 4.4: Problem definitions for validation study 


Problem 

number 

Material 

code 

Aspect 

m 

Thickness 

ratio(a//i) 

Boundary 

conditions 

laminae lay-up 

sequence 

P13 

M3 

1.0 

100 

SSSS 

isotropic 

P14 

M3 

1.0 

100 





1.0 




P16 

M2 

1.0 


SSCS 

[0/90]5 




10 


[0/90]5 

P18 

M2 

1.0 

5 


[0/90]5 

P19 

M2 

1.0 

5 

SSCS 

[0/90]5 

P20 

M2 

1.0 

5 

cscs 

[0/90]5 


The comparative study of the analytical results from the present investigation with 
the results available in the published literature is presented in Table 4.5 . 

From Table 4.5 it can be observed that the difference in the results obtained from 
the present study and that from the published literature is less than 2.5% in all the 
cases. Hence we can conclude that the present formulation gives results which are in 
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good agreement and validates the present formulation. 


Table 4.5: Comparison of results 


Problem 

Laminate theory 

Reference 

Present 

number 

[ref. no.] 


Ny 

Nx 

Ny 

P13 

CLT [16] 

4.0 

- 

3.997 

- 


HSDT 

- 


- 

> 

P14 

CLT [16] 

10.15 

_ 

10.032 

- 


HSDT 

- 


- 

- 

P15 


H 










P16 

CLT [9] 

- 

59.288^ 

- 

60.563" 


HSDT [9] 

- 

m 

- 

33.969 

P17 

CLT [9] 

- 

89.770^ 

- 

91.015" 


HSDT [9] 

- 


- 

35.958 

P18 

CLT [9] 

- 

35.232^ 

- 

36.198" 


HSDT [9] 

- 



12.603 

P19 

CLT [9] 

- 

59.288^ 

- 

60.563" 


HSDT [9] 

- 

12.800^ 

- 

12.969 

P20 

CLT [9] 

- 

89.770^ 

- 

91.015" 


HSDT [9] 

- 

13.659'^ 

! 

13.426 


^ Results obtained with the exact solution developed in the cor- 
responding reference. 

Results using the finite element solution in the corresponding 
reference. 

•= Results obtained by subpressing Ws over the entire domain of 
the plate. 

Nx/Ny are the dimensionless initial buckling values where the 
nondimensionalizing factor is as given in the corresponding 
references. 
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4.6 PARAMETRIC STUDY OF COMPOSITE LAM- 
INATES SUBJECTED TO UNIAXIAL IN- 
PLANE LOADINGS 


In the following section the effects of transverse shear deformation, plate aspect 
ratio (a/6), length/thickness ratio (a//i) and fibre orientation angle (9) on the ini- 
tial buckling loads of laminated (8 layer antisymmetric, [i 0 ] 4 ) composite plates is 
studied (using HSDT) keeping the same material properties (Material- M2). 

The loads considered are uniaxial in-plane compressive and negative shear buckling 
load whereas the geometric boundary condition is taken as simply supported on all 
four sides. 

Results are also obtained using CLT for thin laminate {a/h = 100). 

The results presented in Figs. 4.1 - 4.3 show that for thin laminate {a/h = 100) the 
fibre orientation angle for the maximum value of the fundamental initial buckling 
load (henceforth the fibre orientation angle for the maximum value of the fundamen- 
tal initial buckling load is denoted as 9cric. and the corresponding initial buckling load 
as the critical initial buckling load) is around 45° for aspect ratios a/b = 1,2, 3.The 
critical initial buckling load decreases very slightly with the increase in aspect ra- 
tio.Purther comparing the graphs for a/h = 100 in Figs. 4.1 - 4.3 (obtained using 
HSDT) with that of the graphs in Fig. 4.4 (obtained using CLT) it can be con- 
cluded that the effect of neglecting transverse shear in the analysis is negligible for 
thin laminates. The initial buckling values for thin laminate using HSDT are slightly 
less than the corresponding values obtained using CLT. 

With increase in the thickness of the laminate (decreasing the length/thickness ra- 
tio) transverse shear effect plays a dominant role in the buckling behaviour of the 
laminate.For a/h = 50 the fibre orientation angles for critical initial buckling load 
are about 45°, 42° and 40° respectively, for aspect ratios a/b = l,2,3.Hence the 
critical fibre orientaion angle for a/h = 50 is lesser than that for a/h = 100 and 
changes slightly with the change in aspect ratio. 
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changes slightly with the change in aspect ratio. 

As the laminate becomes thicker there is a shift in the critical fibre orientation 
angle towards lesser values. Also this shift is more as the aspect ratio increases. For 
thick laminate {a/h = 20, 10, 5) the effect of fibre orientation angle on the value of 
critical buckling load decreses with the increase in aspect ratio. Further, it can be 
observed that as the length/ thickness ratio decreases (the laminate becomes thicker) 
the initial buckling load changes very little with the change in fibre orientation angle 
(for all aspect ratios). 

The decrease in the values of initial buckling loads (for any orientation of the fibre) is 
more pronounced for thick laminates. This brings out the fact that for thick laminates 
one should consider transverse shear effects. 

From Fig. 4.5 it can be observed that for shear buckling load (considered negative 
here) the critical fibre orientation angle changes slightly with the change in aspect 
ratio. This critical fibre orientation angle is around 45° for thick laminates {a/h = 
10 ). 

The critical fibre orientation angles and associated critical initial buckling loads 
(uniaxial in-plane compressive, Nx) for different aspect ratios and length/thickness 
ratio are presented in Table 4.6 . 

In Figs. 4.6 - 4.27 the mode shapes of the buckled plates are shown for various 
aspect ratios and length/thickness ratios at different fibre orientation angles.The 
mode shapes are non-dimensionalized as to = ^ where w = tOb-t-Wj.The mode shapes 
are found to be either symmetric or antisymmetric in nature.lt can be observed from 
the figures that there is a change in the mode shapes as the fibre orientation angle 
changes. 
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Table 4.6: Effect of aspect ratio and length/thickness ratio upon the critical 
initial buckling load of antisymmetrically laminated angle ply plates, [±0]4 


Aspect 

Thickness 

Critical fibre 

Critical initial 

(a/6) 

ratio(a//i) 

orientation 

angle(0crtc.) 

buckling load(iV 3 . = N^b^/Erh^) 



HSDT 

CLT 

HSDT 

CLT 

1 

100 

45 

45 

67.04 

67.54 

1 

50 

45 

45 

65.59 

67.54 

1 

10 

42 

45 

40.12 

67.54 

1 

5 

15 

45 

14.66 

67.54 

2 

100 

45 

45 

65.56 

67.52 

2 

50 

42 

45 

61.02 

67.52 

2 

20 

35 

45 

36.65 

67.52 

2 

10 

16 

45 

14.63 

67.52 

2 

5 

15 

45 

4.88 

67.52 

3 

100 

45 

45 

63.195 

67.48 

3 



45 

52.02 

67.48 

3 

20 

27 

45 

22.13 

67.48 

3 

10 

10 

45 

7.64 

67.48 

3 

5 

0 

45 

2.65 

67.48 
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0 10 20 30 40 50 60 70 80 90 

orientation angle 6 (in degrees) 

Figure 4.1: Non-dimensionalized initial buckling load versus fibre orienta- 
tion ( 0 ), (Material- M2). 



orientation angle 0 (in degrees) 

Figure 4.2: Non-dimensionalized initial buckling load versus fibre orienta- 
tion ( 0 ), (Material- M2). 


g] 





tion ( 9 ), (Material- M2). 



Figure 4.4: Non-dimensionalized initial buckling load versus fibre orienta- 
tion ( 9 ) using CLT (Material- M2). 
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Figure 4.5: Non-dimensionalized initial shear(-ve) buckling load versus fi- 
bre orientation (0), (Material- M2). 



Figure 4.6: Mode shape of antisymmetric angle ply laminate [±^]4 for a/b - 
1, a/h = 100, 0 = 30 under application of iV*. 
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nondim. contour lines for w=Wh+Ws — i — 

200 

100 

0 • * 

w=Wt,+Ws(nondim ) -100 



Figure 4.7: Mode shape of antisymmetric angle ply laminate [±0]4 for a/b - 


l,a/h = 100 , 9 = 70 under application of N^. 


w=Wtj+Wg(nondim ) 


nondim. contour lines for w=Wb'*'^s 


5 

0 - “ 



Figure 4.8: Mode shape of antisymmetric angle ply laminate for a/b 
l,a/h = 5,e = 10 under application of iV^. 
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Figure 4.11: Mode shape of antisymmetric angle ply laminate [±d]i for 
a/b = 2,a/h = 100, d = 75 under application of N^. 



Figure 4.12: Mode shape of antisymmetric angle ply laminate [±0]i for 
a/h = 2, a/h = 20,6 = 20 under application of Nx- 
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nondim contour lines for w=W}j+Wg 
20 


w=Wij+Ws(nondim ) 



Figure 4.13; Mode shape of antisymmetric angle ply laminate [±9]^ for 
a/h = 2,a/h = 2Q,d = 75 under application of N^. 



Figure 4.14: Mode shape of antisymmetric angle ply laminate [±^]4 for 
a/6 = 2, a/h = 5,9 = 10 under application of Nx- 
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Figure 4.191 Mode shape of antisymmetric angle ply laminate [±9]. for 
ajh = 3, ajh = 20,0 = 50 under application of iV„. 



» M... .. — “ 

„/i, = 3 . o/h = 5, e = 0 under application of W.. 
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w=Wtj+Ws(nondim ) 


nondim. contour lines for wsw^^+Wg — 
2 — 
1 — 
0 » 
-1 

-2 - • 



Figure 4 . 21 : Mode shape of antisymmetric angle ply laminate [± 0]4 for 
a/h = 3 , a/h = 5, 0 = 40 under application of 



Figure 4.22: Mode shape of antisymmetric angle ply laminate [±^]4 for 
o/6 = 1, a/h = 10, 0 = 30 under application of 
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Figure 4.23: Mode shape of antisymmetric amgle ply laminate [±^]4 for 
a/6 = 1, a/h = 10, 0 = 75 under application of~Nxy. 



Figure 4.24: Mode shape of antisymmetric angle ply laminate [±0]4 for 
a/b = 2, a/h =10,9 = 30 under application of -N^^y. 
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Figure 4.25: Mode shape of antisymmetric angle ply laminate [±0]4 for 
a/b = 2, a/h = 10,9 = 75 under application of -Nxy. 
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Figure 4.26: Mode shape of antisymmetric angle ply laminate [±^]4 for 
a/b = 3, a/h = 10,9 = 30 under application of-N^r^ 

^ ,, .:. .‘'.'..'la '4f4."n 

51 143457 ^^^ 
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Figure 4.27: Mode shape of antisymmetric angle ply laminate [±0]4 for 
a/h — 3, a/h = 10,9 — 75 under application of N^y. 

4.7 INTERACTION CURVES FOR BIAXIAL IN- 
PLANE LOADINGS 

In the present investigation stability boundaries are obtained for two cases of biaxial 
loading. 

• In-plane compressive loads in both X and F directions (i.e. Nx and Ny). 

• In-plane compressive loads in X direction and in-plane (both negative and 
positive) shear (i.e.Nx and i:Nxy). 

Computations are carried out for laminates simply supported on all the four edges 
(SSSS). 

It can be noted here that in this section Nx, Ny and Nxy all denotes nondimen- 
sionalized in-plane initial buckling loads, nondimensionalized by a factor {h’^IErh?). 
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4-7.1 BIAXIAL IN-PLANE LOADINGS iV^ AND Ny 

First the results from the present study are compared with standard results available 
in the literature [15]. It is found that the results from the present study are in good 
agreement with the results cited in the reference. These are shown in Fig 4.28 (for 
a/b = 1) and Fig 4.29 (for a/b = 2) .It can be observed, that, for rectangular 
laminates (in Fig. 4.29) there is a unique point (at about ^ = 0.2667) at which for 
both the laminates ([45/90/ -- 45/90]2, and [±45]4,) the buckling load is same.For 
the combinations of loads where is less than this critical value [± 45 ] 4 s laminate 
buckles at a lower value of Ny for a given Nx and for the combinations of loads where 
Ny is less than this critical value [± 45 ) 44 , laminate buckles at a greater value of Nx 
for the same amount of A^j,.For square laminate no such unique point is observed. 

Next the effect of neglecting transverse shear is studied for square laminated plates.lt 
is found, that, for thin laminates the effect is negligible.lt is observed, that, for 
a/b = 1 (Figs. 4.30 and 4.31) the values of initial buckling load obtained using 
CLT are somewhat more than that obtained using HSDT, but this difference is very 
small and can be neglected for all practical purposes.However, for thick laminates 
the effect is not the same. Here the initial buckling loads obtained using CLT are 
far greater than that obtained using HSDT (Figs. 4.32 and 4.33). Hence if one uses 
CLT to estimate the initial buckling load for thick laminates, it will be grossly 
erroneous.So in this case one has to use HSDT to obtain accurate values of initial 
buckling load. 

In the following section the stability envelops are shown for both cross ply and angle 
ply symmetric laminates [± 0 ] 2 « for plate aspect ratios 1 and 2 and length/thickness 
ratios 100 and 10 for both square and rectangular laminates.The analysis is done 
using simple HSDT for all the cases. 

For a/6 = 1 and a/h = 100 (Fig. 4.34) the following observations can be made. 

• The stability envelops in general are either linear or piecewise linear. 

• ^ = 45° laminate gives the highest value of {Nx + Ny) for any combination of 
the loads and for cross ply laminate {Nx ± Ny) is the minimum. 
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• There is a sudden change in slope of the stability envelops for laminates with 
0 = 30 and Q = 60°. For Q = 30° laminate, this slope change occurs at small 
value of and for B = 60° laminate slope change occurs at high value of 
A^aj.For 0 = 45 laminate and for cross ply laminate there is no change in the 
slopes of their stability envelops. 

• For Q = 30° and B = 60° laminates the portion where the slopes of the two 
envelops are the same, the values of (iV^ + Ny) are also same. 

• From the stability envelops shown in Fig. 4.34 no such unique point, as men- 
tioned earlier (Fig. 4.29) can be observed. Instead here it can be observed that 
for B = 30° and B = 60° laminates there axe two points (approximately at 
^ = 0.175 and ^ = 0.175) between which both the laminates have the same 
interaction curve For the combinations of loads where is less than the crit- 
ical value ( ^ = 0.175) 0 = 60° laminate buckles at a greater value of Ny 
for the same amount of Nx and for the combinations of loads where Ny is less 
than the critical value ( ^ = 0.175) ^ = 30° laminate buckles at a greater 
value of Nx for the same amount of iVj,. Stability envelops of both the cross ply 
laminate and 6 = 45° laminate donot intersect any any other curves. 

For a/b = 1 and a/h = 10 (Fig. 4.35) the following observations can be made.There 
is a change in the nature of the stability envelops for thick square laminate from 
that of the thin laminate though the basic shapes of the stability envelops remains 
the same (linear or piecewise linear curves). 

• Except for cross ply laminate (linear stability envelop) the sability envelops 
for the angle ply laminates are all piecewise linear. 

9 d z= 45° laminate gives the highest value of {Nx + Ny) for almost all combi- 
nations of biaxial load except for very high and very low values of Nx and for 
cross ply laminate {Nx + Ny) is the minimum for any combination of biax- 
ial loading. For very high values of Nxi 0 = 30° laminate gives the maximum 
value of {Nx + Ny) and for very low values of Nx, 0 = 60° laminate gives the 
maximum value of {Nx + Ny). 
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. Change in slope for both 9 = 30" and 9 = 60" laminates occur when both N, 
and N, are moderate.For « = 45" laminate change in slope occurs both at very 
high and very low values of N,.ln all the cases a change in slope is aasociated 
with a corresponding change in mode shape which is evident from Figs. 4.38 
- 4.45. 


• The portion for which + Ny) is same for both e = 30° and 0 = 60° 
laminates is less as compared to thin laminate. 

• Here again, observations, similar to that for thin laminates regarding inter- 
section points of the stability curves of different laminates, can be made for 
9 = 30° and 9 = 60° laminates.Unlike for thin laminate here 9 = 45° laminate 
cuts the 9 = 60° laminate at very high value of Ny and the 9 = 30° laminate 
at very high value of N^- 

The stability boundaries for rectangular laminated plates {a/b — 2) are quite differ- 
ent from that of square laminates. 

For a/b = 2 and a/h = 100 (Fig. 4.36) the following observations can be made. 

• Except for 9 — 30° laminate (linear stability envelop) the sability envelops of 
other angle ply laminates and cross ply laminate are all piecewise linear. 

• For different combinations of Nx and Ny, either 6 = 45° laminate oi 9 = 60° 
laminate gives the maximum value of {Nx + Ny). 

• For cross ply laminate, the slope of the stability envelop changes at a high 
value of Ny, for 9 = 45° laminate change in slope occurs at a high value of Nx 
whereas for 9 = 60° laminate change in slope occurs twice.One is when both 
Nx and Ny are high and another is when only Nx is high and Ny is moderate. 

• Here four points can be observed at which different laminates (two laminates 
at each points) gives the same value of (iVt + iVj,).That is at these points the 
biaxial load carrying capacities of two different laminates are the same. 

For o/6 = 2 and a/h = 10 (Fig. 4.37) the following observations can be made. 
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• Here the stability envelops are piecewise linear but for certain portions of the 
envelops, slope changes very frequently. 

• The curves arc almost perpendicular to the V axis (denoting Nyb'^/Erh^) for 
one portion of them and are almost perpendicular (except for 9 = 45° laminate) 
to the X axis (denoting /Exh^) for another portion.This manifests to the 
fact that at these portions of the curves, for a small decrease in the value of 
Nxi the capacity of the laminates to carry Ny increases rapidly and vice-versa. 

• Change in slopes of all the envelops occur at high values of A(c. Similar to a 
square laminate here also a change in slope is is aasociated with a correspond- 
ing change in mode shape which can be observed from Figs. 4.46 - 4.57. 

• As that for thin rectangular laminates, similar kind of intersection points 
can also be observed for thick rectangular laminates.Here the number of such 
unique points are six. 
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Figure 4.28: Comparative study of Stability envelop for biaxial in-plane 
compressive loads and Ny, (Material- M4), (6 = 254mm, h = 2.112m7n), 
from [15]. 



Figure 4.29: Comparative study of Stability envelop for biaxial in-plane 
compressive loads Nx and Ny, (Material- M4), (i> = 254, h = 2.112mm), 
from [15]. 
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a/h=100 (using HSDT 
a/h=100 (using W 5=0 


a/b=1 

For[ 30 /- 30 ] 2 s 




Figure 4.30: Effect of transverse shear on stability envelop for biaxial in 
plane compressive loads Nx and Ny, (Material- M2). 
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Figure 4.31: Effect of transverse shear on stability envelop for biaxial in 
plcine compressive loads Nx and Ny, (Material- M2). 
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Figure 4.32: Effect of transverse shear on stability envelop for biaxial in 
plane compressive loads and Ny, (Material- M2). 


a/h=10 (using HSDT) 
a/h=10 (using Ws=0) 


For [60/-60]2 s 




Figure 4.33: Effect of transverse shear on stability envelop for biaxial in 
plane compressive loads Nx and Ny, (Material- M2). 







Figure 4.34: Stability envelop for biaxial in-plane compressive loads and 
Ny, (Material- M2). 



Ny/E,h® 

Figure 4.35: Stability envelop for biaxial in-plane compressive loads and 
Ny, (Material- M2). 
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Figure 4.36: Stability envelop for biaxial in-plane compressive loads and 
Ny, (Material- M2). 



Figure 4.37: Stability envelop for biaxial in-plane compressive loads No; and 
Ny, (Material- M2). 
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Figure 4.38: Mode shape of symmetric cross ply laminate [90/0]2s for a/b 
l,a/h = 10 at ^ = 0.6. 
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Figure 4.39: Mode shape of symmetric angle ply laminate [±30]2s for a/b 
l,a/h = 10at ^ = 0.2. 
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Figure 4.40: Mode shape of symmetric angle ply laminate [±30]2s for ajh - 
l,a//i = 10at ^ = 0.4. 



Figure 4.41: Mode shape of symmetric angle ply laminate [±45]2. for ajh 
l,a//i = 10 at ^ = 0.1. 
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Figure 4.42: Mode shape of symmetric angle ply laminate [±45]2s for a/h 
l,a//i = 10 at ^ = 0.8. 
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Figure 4.43: Mode shape of symmetric angle ply laminate [±45]2. for ajh 
1, ajh = 10 at ^ = 0.1. 
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Figure 4.48: Mode shape of symmetric cross ply laminate [90/0]2s for a/6 — 
2,a//i = 10 at ^=0.2. 



Figure 4.49: Mode shape of symmetric angle ply laminate [±30]2s for a/6 - 
2, ajh = 10 at ^ = 0.8. 
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Figure 4.50: Mode shape of symmetric angle ply laminate [±30]2s for ajb = 
2,a//i = 10 at ^ =0.2. 



Figure 4.51: Mode shape of symmetric angle ply laminate [iSOja* for o/6 - 
2,a//i=10at ^ = 0.1. 
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Figure 4.54: Mode shape of symmetric angle ply leiminate [±45]2s for o,jh = 
2,a//i= 10 at ^ =0.2. 


Figure 4.55: Mode shape of symmetric angle ply laminate [±60]2s for a/h = 
2,o//i = 10 at ^ = 1.0. 
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4.7.2 BIAXIAL IN-PLANE LOADINGS AND ±N^y 

First the results from the present study are compared with standard results available 
in the literature [15]. It is found that the results of the present study are in good 
agreement with the results in the reference.These are shown in Fig 4.58 (for a/b = l) 
and Fig 4.59 (for a/b — 2).It can be observed, that, for rectangular laminates (in Fig. 
4.59) there is a unique point (at about = 0.148) at which both the laminates 
([45/90/ - 45/90] 2 s and [±45]4 s) carry the same amount of axial compressive and 
shear loads, and NxyPdso when only positive shear load is applied both the 
laminates carry the same amount of loads whereas the same is not true for negative 
shear load. For square laminate no such obseravations can be made. 

In the following section stability envelops are shown for both cross ply and angle ply 
symmetric laminates [±^] 2 s with plate aspect ratios 1 and 2 and length/thickness 
ratios 100 and 10 for both square and rectangular laminates.The analysis is done 
using simple HSDT for all the cases. 

For a/b = 1 and a/h = 100 (Fig. 4.60) the following observations can be made. 

• Stability envelops are not symmetric about N^y = 0 line except for the cross 
ply laminate.lt is observed that the value of {Nx + \Nxy\) for initial buckling 
load is more for negative shear as compared to positive shear. 

• For all combinations of biaxial loading cross ply laminate gives the lowest value 
of {Nx + \Nxy\) and 6 = 45° laminate gives the highest value of [Nx + liVa;j,|). 

• Stability envelops of 6 = 30° and 9 = 60° laminates are comparable though 
9 = 30° laminate gives slightly higher value of {Nx + \Nxy\). 

• There is no such unique point for biaxial loading at which two or more lami- 
nates carry the same amount of axial compressive and shear loads.When only 
shear load is applied 9 = 30° and 9 = 60° laminates carry the same amount of 
load. 

For a/5 = 1 and a//i = 10 (Fig. 4.61) the following observations can be made. 
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• Stability envelops are not symmetric about Nxy = 0 line except for the cross 
ply laminate.lt is observed that the value of {N^ + for initial buckling 
load is more for negative shear as compared to positive shear. 

• For positive shear loading the stability envelops of 0 = 30° and 9 = 45° lam- 
inates are comparable with each other while 9 = 60° and cross ply laminates 
are comparable with each other. 

• For negative shear loading the stability envelops of 0 = 30° and 9 = 45° 
laminates are comparable with each other but {Nx + \Nxy\) is more for 9 = 60° 
laminate than for cross ply laminate. 

• Here the stability envelops of different laminates intersect each others at some 
unique points. At these points different laminates carry the same amount of 
biaxial loads, Nx and Nxy.ln other words it can be said that there are some 
unique combinations of in-plane compressive and shear loadings that can be 
carried by different laminates. 

For a/6 = 2 and a/h = 100 (Fig. 4.62) the following observations can be made. 

• Stability envelops are not symmetric about Nxy = 0 line except for the cross 
ply laminate.lt is observed that the value of {Nx + jAij,!) for initial buckling 
load is more for negative shear as compared to positive shear. 

• For all combinations of biaxial loading except for very high values of Nxy (both 
positive and negative), 9 — 45° laminate gives the highest value of (A(e+|-^sj/|)- 

• The range of shear loading (defined as: positive shear - negative shear) which 
a laminate can withstand is minimum for 9 = 45° laminate. 

• Similar kind of obseravtions that were made for square thin laminates can also 
be made for rectangular thin laminates regarding intersection poins of stability 
envelops. 

For a/b = 2 and a/h = 10 (Fig. 4.62) the following observations can be made. 
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• Stability envelops are not symmetric about N^y = 0 line except for the cross 
ply laminate.lt is observed that the value of (iV^ + \Nxy\) for initial buckling 
load is more for negative shear as compared to positive shear. 

• For low values of shear load (both positive and negative), the difference in 
{Nx + |iVxy|) is substantial for different laminates whereas for high values of 
positive shear load this difference is very nominal. For very high values of neg- 
ative shear load the above mentioned difference is medium. 

• For negative shear load cross ply laminate gives the minimum value of (Nx + 
\Nxy\) except for very low values of negative shear load. 

• For high values of Nx, 0 = 30° laminate gives maximum value of (Nx + 

• Here also one can observe some unique points at which different laminates 
carry the same amount of biaxial loads, Nx and Nxy 


* In the above discussion denotes the absolute value of nondimensionalized 

in-plane shear buckling load. 
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Figure 4.58: Comparative study of Stability envelop for biaxial in-plane 
compressive and shear loads and N^y, (Material- M4), (b = 254mm,/i = 
2.112mm), from [15]. 



compressive and shear loads Nx and Nxy, (Material- M4), (i> — 254mm, h — 
2.112mm), from [15]. 
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Figure 4.60: Stability envelop for biaxial in-plane compressive and shear 
loads Nx and N^y, (Material- M2). 
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Figure 4.61: Stability envelop for biaxial in-plane compressive and shear 
loads Nx and Nxy, (Material- M2). 
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Figure 4.62: Stability envelop for biaxial in-plane compressive and shear 
loads Nx and Nxy, (Material- M2). 



Figure 4.63; Stability envelop for biaxial in-plane compressive and shear 
loads Nx and Nxy, (Material- M2). 




Chapter 5 


CONCLUSIONS AND FUTURE 
SCOPE 

5.1 CONCLUSIONS 

In the present work initial buckling behaviour of composite laminated plates was 
carried out with both uniaxial and biaxial loads (both compressive and shear) in- 
cluding transverse shear effect.The response has been obtained for different plate 
aspect ratios, fibre orientation angles and length/ thickness ratios.The analysis is 
performed using finite element method and the eigenvalues (which gives the initial 
buckling load) and the corresponding eigenvectors (which gives the mode shape for 
the corresponding buckling load) are obtained from the finite element solution us- 
ing standard NAG routines. The following main conclusions can be drawn from the 
present investigation. 

1. For thin laminates the effect of transverse shear is negligible and CLT is suf- 
ficient for predicting buckling loads. However as the length/thickness ratio increases 
transverse shear plays a dominant role and it has to be duly accounted for accu- 
rate prediction of critical loads..For moderately thick laminates there is substantial 
change in the values obtained by CLT and HSDT.For accurate analysis HSDT must 

be used for thick/moderately thick laminate. 

2. The effect of lamina orientation angle of antisymmetric laminate on the initial 
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buckling load decreases with the increase in length /thickness ratio. 

3. For thick laminates, the effect of lamina orientation angle of antisymmetric 
laminate on the initial buckling load decreases as the aspect ratio increases. 

4. The critical fibre orentation angle reduces as the length/thickness ratio of the 
laminate increases. 

5. For thick laminate the critical fibre orentation angle becomes smaller as the 
aspect ratio of the laminate increases. 

6. For shear buckling load, as the aspect ratio increases the effect of fibre orien- 
tation angle on the value of initial buckling load decreases substantially. 

7. The critical value of the fibre orientation angle remains almost same for shear 
buckling load with change in aspect ratio. 

8. Interaction curves for Nx and Ny are either linear or piecewise linear. 

9. There are some combinations of Nx and Ny for which different laminates give 
same value of critical load. 

10. For combined shear and compressive load {Nxy and Nx) the stability envelops 
are not symmetric about the the line Nxy = 0 except for cross ply laminate.For same 
Nx lamiante can withstand more negative shear loading than positive one. 

11. For some combinations of combined compressive and shear loadings, the 
stability envelops for different laminates are culsered and for some other portions 
they are widely separated. 

5.2 SCOPE FOR FUTURE WORK 

Some of the areas which need to be investigated further are mentioned below. 

1. The analysis for laminates with cut-outs of different geometries (practical 

laminates). 

2. Response of laminates with plydrops need to be investigated. 

3. Stability characteristics of thick laminates under the application of partial 
edge loading and concentrated loading can also be looked into. 

4. Post-buckling analysis of laminates using non-linear strain-displacement velar 

tions. 
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5. Stability envelops of thick laminates with different combinations of ther- 
mal/hygroscopic and static in-plane loading. 
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APPENDIX B 


The element shape functions of a four noded rectangular continuous finite element 


in local co-ordinates are: 
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The aforementioned shape functions can be expressed in natural co-ordinates as: 

- 2) (rj - If {-ri - 2) 

^2 = + Ifu - 2)(r; - - 2) 
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The details of [Sq] and [Bg] matrices are as shown: 
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In the above matrices, for columns 1 and 2, i = 1 — 12 (shape functions cor- 
responding to Mo and vq) and for columns 3 and 4, ^ = 1 — 16 (shape functions 
corresponding to and Wy). 
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